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Abstract 

We use a new version of the reversible Parking Lot Model to study the compaction of vibrated 
polydisperse media. The particle sizes are distributed according to a truncated power law. We 
introduce a self-consistent desorption mechanism with a hierarchical initialization of the system. 
In this way, we approach densities close to unity. The final density depends on the polydispersity 
of the system as well as on the initialization and will reach a maximum value for a certain exponent 
in the power law. 



I. INTRODUCTION 

The vibratory compaction of granular ma- 
terials has long been of importance in tech- 
nological applications like high performance 
concrete or ceramics which have to with- 
stand extreme stress. Extensive experimen- 
tal studies have been conducted on monodis- 
perse systems. Knight et al. Q] investi- 
gated density relaxation of a column contain- 
ing monodisperse spherical beads subject to 
a long sequence of taps while Nowak et. al P] 
studied density fluctuations in vibrated gran- 
ular materials. 

Several models were introduced to de- 
scribe the dynamics of a granular system un- 
der compaction, including a lattice model [3J] 
and the Tetris model special case of 
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a frustrated lattice gas. The study of poly- 
disperse systems reaches back to the ancient 
Greeks where Apollonius of Perga studied the 
problem know as 'Apollonian packing' P, f|. 
This problem deals with the question how 
to tile the space with circles by iteratively 
placing, between every three circles, a circle 
tangentially touching all three. The Apollo- 
nian packing is a special case of the so called 
'space filling bearings' (SFB) in which a plane 
is tiled with circles touching one another such 
■;hat the whole area is covered with circles 
3, 0, ^, 10 1 . These space filling bearings fill 



space with a particle size distribution given 
by a truncated power law. Different stud- 
ies of polydisperse packings where conducted 
by Aste (llj , Dodds and Weitz 



liantov et al. 113, 1 



12 and Bril- 



The aim of our work is to study the time state with a final density of Pj am ~ 0.7475. 

evolution of density in polydisperse systems JZ •, , -o AT . , , . 

J Krapivsky and Ben-Maim extended this 

under vibratory compaction. The size dis- , , , -ut t-> 1 • t ± u j 1 

J model to the reversible Parking Lot Model 
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tribution of the particles obeys a truncated L J TT . . , .,, , , 

j 1 olj . Here particles will additionally desorb 

power law (different than the one used by ., , , . . . /r . m, , , , 

v with a desorption rate k [r igure UJ) and the 



Brilliantov) . We modified the one dimen- 



sional reversible Parking Lot Model 
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density of the system's equilibrium state can 

be different from Pj am . This equilibrium state 

171 where identical particles adsorb on an in- . , , ,, , . , , , , , 

— 1 ^ is determined by the ratio between the ad- 

terval with an adsorption rate k + and desorb . . , , . , , i+ n - 

sorption and the desorption rates k = k^/k . 

with a rate k~ , such that the system is hierar- 



chically initialized and the equilibrium state 
with the highest density is reached through 
self-consistent desorption. 

II. THE MODEL 

We first explain the classical Parking Lot 
Model 3 and its reversible variant Il5|]. Af- 
terwards we introduce our generalization to 
polydisperse systems and the self consistent 
desorption probability. 
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FIG. 1: The reversible Parking Lot Model with 
an interval length I and particles of equal size r. 
Particles adsorb with a rate k + and desorb with 
a rate k~ . 



For the final density p ss the following leading 
A. The Parking Lot Model behavior in the two limiting cases was found: 

In the classical Parking Lot Model fPLM), 

!k for k <^ 1 

(1) 
1 — ln(k) for k ^> 1 

places along an interval of length I. The ad- 
sorption fails if the chosen place is partially In the limit of k — > oo (k + = l,k~ — > 0), also 
occupied by a previously adsorbed particle, called the desorption controlled limit, the in- 
Because of the irreversible adsorption mecha- terval is completely filled with particles and 
nism the system will reach a so called jammed p ss = 1 when time goes to infinity. 



B. Generalization to a polydisperse 
Model 

Granular materials that reach highest den- 
sities consist of particles with different sizes. 
Such media are called polydisperse. We will 
simulate such media using a polydisperse 
PLM. 

The behavior of a polydisperse PLM is 
governed by the size distribution of the parti- 
cles, the amount of particles available in the 
system, the initialization process and the de- 
pendence of the desorption rate on the par- 
ticle size. Figure E] shows the setup for the 
polydisperse Parking Lot Model. 



hi * | p(r) N^~~^/ 
I Reservoir 

FIG. 2: The self consistent reversible polydis- 
perse Parking Lot Model with an interval length 
Z, particles of different size r and the reservoir. 
Particles adsorb with a rate k + and desorb with 
a conditional probability p(r) which also de- 
pends on the distribution of holes in the interval. 

Unlike to the case of monodisperse sys- 
tems we introduce a reservoir because it turns 
out to be necessary to restrict the number 
of particles of each size in order to impose a 
given size distribution. For a system with K 
particles of which K 1 are adsorbed on the in- 
terval, this reservoir will contain the remain- 



ing K 2 — K — K\ particles. 

We restrict the size of the particles to the 
interval R = [r min , r max ] which we will re- 
fer to as the 'range of particle sizes'. For 
all simulations discussed here, r max is set to 
unity and the size distribution of all particles 
is given by a truncated power law: 

P(r)=a-r-\ r eR (2) 

Where a is a normalization factor chosen 
such that fl' max P(r)dr = I. We allow the 
exponent b to take values in the range of 
1.306 < b < 1.802, the same range as found 
for two dimensional space filling bearings . 
We also apply periodic boundary conditions 
to the system. Its density p is given by. 

P = X>i/J (3) 

i=i 

where K x is the total number of adsorbed 
particles and is the length of the zth ad- 
sorbed particle. 

As in the classical PLM, the particles at- 
tempt to be adsorbed with the adsorption 
rate k + . This rate is defined as the number of 
particles attempting to be adsorbed per time 
unit. A particle of given size r, randomly cho- 
sen among the adsorbed particles, will des- 
orb with a conditional desorption probability 
p(r). This conditional probability can be cal- 
culated through: 

p{r) = j:\h-r) /I (4) 
i=i 



Here K\ is the total number of holes which 
is the same as the number of adsorbed par- 
ticles when using periodic boundary condi- 
tions. The variable hi is the length of a sin- 
gle hole, which can be zero in the case of two 
touching particles. The primed sum only con- 
siders terms with hi > r. Thus p(r) will be 
zero if the particle size r is larger than any 
available hole. We call p(r) the 'self consis- 
tent desorption probability' because of its de- 
pendence on the hole distribution. It changes 
in time as the configuration in the interval 
changes. 

This approach is justified because in an ex- 
periment, a particle can only leave its place if 
there is a large enough hole where it can move 
to. The holes of the system are stored in a 
list H which is updated after each adsorption 
and desorption. 

In the monodisperse Parking Lot Model 
a constant desorption rate removes particles 
regardless of their local environment. There- 
fore, the system ultimately reaches an equi- 
librium state independent of the initial con- 
dition. This is not the case in a polydisperse 
system where the self consistent desorption is 
implemented. If we would use a constant des- 
orption rate in the polydisperse system, the 
result would be a final density lower than the 
Pjam obtained in the monodisperse model. 
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1. Dynamics 

Now we will explain the dynamics of the 
system. The word 'random' will denote an 
equally distributed random number genera- 
tor unless stated otherwise. 

First the particles available to the system 
are generated. For this a total number of K 
particles is put into an initially empty reser- 
voir (list). The size of each particle is sam- 
pled randomly from a power law distribution 
Eq. (J2J). The total number K of particles is 
given by J2k=o r k — I- Thus the generation of 
particles is stopped as soon as the sum of the 
lengths of all particles is equal to I. We call 
the size distribution of particles in the reser- 
voir at this time (t = 0) the 't distribution'. 

Next an initial distribution of particles on 
the interval is generated. We will call this 
the 'initialization of the interval'. This is nec- 
essary because for an almost empty interval 
p(r) is close to unity for all particles and any 
adsorbed particle would most likely desorb 
again. A further justification for this initial- 
ization is given in section IIII Al 

In order to reach highest densities a hi- 
erarchical initialization is used. This idea 
was inspired during a visit to W. Losert's 
|l8| laboratory at the University of Mary- 
land. The reservoir is sorted by size such 
that r\ > r 2 > ■ • • > tk- Then each par- 
ticle, starting with the biggest particle (ri), 



is given / trials to adsorb. During the ini- 
tialization, the desorption probability is set 
to zero. This resembles the filling of a recip- 
ient in an experiment, where the large par- 
ticles are loaded first and the system is not 
vibrated. 

For each adsorption trial, a random point 
in the interval is chosen. If the chosen place is 
not even partially occupied the particle will 
adsorb. If the particle did not adsorb after 
its Jth trial it is left in the reservoir and one 
continues with the next smaller one. The ini- 
tialization is finished when each particle had 
up to I trials for adsorption. The density re- 
sulting from this initialization is called Pmit- 

Now the adsorption/desorption mecha- 
nism is activated. With a rate k + , a parti- 
cle is randomly picked from the reservoir and 
adsorbed if possible on a random position. 
If the adsorption is not possible, the parti- 
cle does not get a second trial but will be 
put back into the reservoir. For desorption 
a particle on the interval is randomly picked 
at each time step. With the probability p(r), 
given by Eq. (j3J) , it is put back into the reser- 
voir. Each time the distribution of particles 
in the interval changes the list of holes H is 
updated accordingly. 

In this model we have two different defini- 
tions of time. The first definition covers the 
initialization of the system. Here one time 
step consists of a single adsorption trial. The 
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second definition considers the time after the 
initialization (the dynamics). Here a single 
time step in the system consists of an ad- 
sorption and a desorption trial. Adsorption 
is attempted with the rate k + while desorp- 
tion is checked every time step. These two 
definitions are of significant difference. Nev- 
ertheless we will put them on the same scale 
in our density plots in order to investigate 
dependencies on different parameters. 

III. MONTE CARLO RESULTS 

In this section we give an overview over 
the performed simulations. We always use 
an adsorption rate k + = 1. Thus, one ad- 
sorption trial is performed during each time 
step. 

The alignment of the particles on the in- 
terval / at a certain time is visualized in a 
spatio-temporal diagram by plotting at each 
time step the system by a height of one pixel 
and a width proportional to I. The density is 
gray scale encoded using white pixels for zero 
density and black pixels for density one. The 
darkness of the pixel depends linearly on the 
density. Arranging these pictures in chrono- 
logical order results in the spatio-temporal di- 
agram as displayed in Figure El which shows 
the time evolution of the system during the 
initialization process. 

Because of the hierarchical initialization 




FIG. 3: Spatio-temporal diagram of the initial- 
ization process. System parameters: I = 300, 
r min = 0.001, r max = 1, k+ = 1, I = 1000, 
b = 1.33, K = 5420. Density at the end of the 
initialization: p ~ 0.985 

the big particles, shown as regions of high 
density, are adsorbed first. As time pro- 
gresses the interval is consecutively filled with 
smaller and smaller particles which settle in 
the remaining gaps until the system reaches 
a high density (desorption is always turned 
off in this case). 

A. Size Distribution 

The hierarchical initialization avoids the 
exclusion of large particles and thus leads to 
an adsorption of almost all particles. In the 
ideal case of a completely filled interval, the to 
distribution and the distribution of adsorbed 
particles for a system in its final state would 
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FIG. 4: Particle size distribution for the self con- 
sistent system (solid line, I = 1000) and a sys- 
tem with a constant desorption rate using no ini- 
tialization (dotted line) . In both cases we display 
the distribution at the end of the simulation. 

be identical. 

In Figure E| we display the size distribu- 
tion of the adsorbed particles at the end of 
the simulation. The solid line shows the dis- 
tribution for a self-consistent desorption rate 
as given by Eq. (J3J at the end of the simu- 
lation which coincides very well with the to 
distribution which is omitted for clarity. All 
simulations show the same overall behavior 
for different values of R and b. In addition we 
show the distribution for a system with con- 
stant desorption rate using no initialization. 
For large particles the distribution deviates 
very strongly from the to distribution which 
justifies our hierarchical initialization. 



0.006 p 
0.005 \ 



-§ 0.004 - 




0.001 0.01 0.1 1 
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FIG. 5: Desorption probability of the self con- 
sistent system ELS £1 function of the particle size 
in a semi-logarithmic plot. For a denser system 
the probability is smaller. The solid line is the 
probability at the end of the initialization pro- 
cess (t ~ 10 7 ). The dotted line shows p(r) at the 
end of the simulation (t ~ 10 9 ). System param- 
eters: I = 5000, r min = 0.001, r max = 1, k + = 1, 
I = 1000, b = 1.33, average K = 89618. Data 
are averaged over 10 runs. 

B. Desorption Probability 

The desorption probability p(r) right af- 
ter the initialization of the interval and at 
the end of the simulation is shown in Figure 
EJ The probability depends on the hole dis- 
tribution and therefore changes in time. The 
solid line in Figure shows the probability at 
the end of the initialization process (t ~ 10 7 ). 
At this time, the transition from the initial- 
ization to the adsorption/desorption process 
takes place. The dotted line shows p(r) at the 
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end of our simulation (t ~ 10 9 ). The prob- 
ability p(r) decreases the denser the system 
gets. 

The probability p(r) will become zero for 
large particles at some time and they will not 
be able to desorb anymore. This leads to a 
stable configuration with highest densities in- 
cluding all the large particles in the interval. 

C. Density evolution 




time 

FIG. 6: Density as a function of time in the self- 
consistent system for different values of /. The 
dotted line is the fitted curve according to Eq. 
(jSJ). Fit parameters: = 1, Apoo = 0.967, 
B = 2.829, r = 19339. Data are averaged over 
10 runs. The system parameters are the same as 
in Figure 03 

The development of the density is shown 
in Figure |H1 Here the density evolution 
for different numbers of trials I are shown. 
The inset enlarges the time interval where 



the transition from the initialization to the 
adsorption/desorption process takes place. 
From this point on, adsorption- and desorp- 
tion take place simultaneously; the density 
rises slower. 

Knight [1] and Nowak measured experi- 
mentally that the time evolution of density in 
monodisperse systems can be asymptotically 
fitted by the expression 



P(t) = Pc 



(5) 



1 + Bln(l + t/r) 

where the parameters p^, Ap^, B and r de- 
pend on the experimental setup. Here p^ is 
the final density while Apoo is the difference 
between this value and the density at t = 0. 
The parameters B and r fit the logarithmic 
behavior of the curve. The same asymp- 
totic behavior was obtained by Krapivsky 
and Ben-Nairn for the reversible parking lot 
model using a variable desorption constant 
Q|. The dotted line in Figure El is the fit us- 
ing Eq. (jSJ). We can see that the density 
after the initialization is very close to the fi- 
nal density. Thus the dominating part in our 
simulations is the initialization process which 
we will investigate closer in the following sec- 
tion. 

D. Discussion 

Using the self consistent desorption, the 
final density mostly depends on the initial- 



ization. Large particles that could not ad- 
sorb during the initialization are, in general, 
excluded from the interval for the rest of the 
simulation. The initialization is character- 
ized by pinit- Three parameters determine 
this density. The number of trials /, the ex- 
ponent in the power law b and the range of 
particle sizes R. 

1. Dependence on I 
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FIG. 7: Dependence of pi n u on the number of 
trials I. Displayed are three curves with differ- 
ent r m i n while r max = 1 is kept constant for all 
three runs. The maximal possible density is de- 
termined by r m i n and reached for very large I. 
System parameters: I = 5000, r max = 1, k + = 1, 
b = 1.33 and K ~ 89000. Data are averaged over 
9 independent runs. 

As already seen in Figure El a larger num- 
ber of trials I increases considerably the final 
density pi n u of the system. In Figure [7| the 
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dependence of pi n u on / is displayed. Start- 
ing from a low density of about 0.55 the den- 
sity of the system increases until it reaches its 
maximum density for this value of r min . Be- 
cause the system is hierarchically initialized 
the final density depends on I. This depen- 
dence can be fitted by: 

Pinit(I) = Pmax(r m in) ~ Ap • I~ f (6) 




10 100 1000 

I (number of trials during initialization) 



10000 



Here p m ax{ r min) is the maximum density the 
system can reach for the chosen value of r min 
in the limit I — > oo and Ap is the differ- 
ence between p maa: (r min ) and the initializa- 
tion density for 1 = 1. The exponent / is a 
fit parameter (/ = 0.4406 for r min = 0.001) 
giving the logarithmic relaxation towards the 
highest density. In Figure H the curve fitted 
with Eq. (JBJ) is the dotted line. 

A better fit can be accomplished when in 
Eq. (jUJ) the last term is multiplied by the 
last term given in Eq. (j3J). The resulting 
expression has the form: 

Ap 



FIG. 8: Simulation data of Figure El fitted by 
Eq. ® (dotted line, p max = 1, Ap = 0.4406 , 
/ = 0.4109) and Eq. © (solid line, p max = 1, 
Ap = 0.4406, B = 0.636, r = 21.3, f n = 0.3163.) 
The system parameters are the same as in Figure 
7 
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(7) 



In Figure |Hl this results in the full line coin- 
ciding with the simulation data. The value 
of /„ is /„ = 0.3163 for r min = 0.001. 

The number of trials I also influences the 
size distribution of the particles. In Figure 
IH1 the size distribution for two different val- 
ues of / is displayed. The distribution of ad- 
sorbed particles at the end of the simulation 



FIG. 9: Comparison of the size distribution at 
the end of the initialization for a different num- 
ber of trials / (dotted line) with the to distribu- 
tion (solid line). System parameters: I = 5000, 



0.001, r r , 



1, k + = 1, b = 1.33 and 



K « 89000. 



coincides better with the to distribution for 



larger /. 

The deviation from the original distribu- 
tion depends on the value of I because a 
larger / gives each particle more trials for 
adsorption. For large J, big particles ad- 
sorb first because of the hierarchical ini- 
tialization and a deviation is only notice- 
able for small particles. Comparing Figure 
with Figure 0] reveals that the deviation 
from the to distribution (large J) after the 
initialization almost vanishes when the ad- 
sorption/desorption mechanism is activated. 

Finally we study which effect I has on the 
self consistent desorption probability p(r). 
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FIG. 10: Desorption probability for different 
number of trials I in a double logarithmic plot, 
right after the initialization. It is zero for the 
largest particles in all three cases. System pa- 



rameters: I = 5000, r« 



0.001, TV 



1. 



k+ = 1, b = 1.33, and K ps 89000. 



In Figure [TUl the desorption probability 
right after the initialization is displayed for 



three different I. The larger I the smaller 
the overall desorption probability. More tri- 
als / during the initialization result in smaller 
holes and a denser system. Hence the overall 
probability will be lower for a larger I. In all 
three runs it is zero for the largest particles. 

2. Dependence on R 
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FIG. 11: Dependence of pi n u on r m i n fitted by 
Eq. © {p m ax(I) = 0.987, Ap = 0.241, e = 
0.600 for I = 1000 and p max {I) = 0.562, Ap = 
0.093, e = 0.519 for 1 = 1). The smallest particle 
size r m i n is varied while r max is kept constant. 
System parameters: I = 5000, r max = 1, k + = 1, 
b = 1.33, K ps 5200 - 90000. Data are averaged 
over 9 runs. 

The maximum density also depends on the 
range of particle sizes R. In Figure [TT] r max 
was kept constant (unity) while r m i n was var- 
ied. The dependence of pi n u on r min is dis- 
played for two different values of /. A larger 
range of particle sizes R results in a higher 
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Pinit- This general behavior is independent of 
the number of trials /. 

The hierarchical initialization, starting 
with the biggest particle, allows the system to 
densify better because the holes get smaller 
and smaller. Would the initialization start 
with the smallest particles first, the system 
would end up in a final state with an ex- 
tremely low density. The smallest size r min 
in the range of particle sizes R determines 
the smallest hole that can be filled and thus 
determines the highest reachable density. 

The dependence on the minimal particle 
size in Figure m can be fitted very well by 
the equation. 



Pinitymin) Pmaxi^I*) Ap ' 



e 

rain 



(8) 



Here p max (I) is the maximal possible density 
for the applied number of trials / which oc- 
curs for r m i n 0. The variable Ap is the 
difference between this maximal density and 
the density the monodisperse system reaches 
for the same applied number of trials /. Thus 
we can determine these values through aver- 
aging over several independent simulations. 
The fit parameter e has a value of e = 0.6 for 
I = 1 and e = 0.519 for I = 1000. 

3. Dependence on b 

The most interesting behavior of the sys- 
tem is revealed when investigating the depen- 




1.3 1.35 1.4 1.45 1.5 1.55 1.6 1.65 1.7 1.75 1.8 1.85 
b (exponent in the power law) 




1.3 1.35 1.4 1.45 1.5 1.55 1.6 1.65 1.7 1.75 1.8 1.85 
b (exponent in the power law) 



FIG. 12: Dependence of Pi n a on the exponent 
b. System parameters: / = 5000, r m i n = 0.001, 
r max = 1, k + = 1, K « 85000 - 422000. Data 
are averaged over 9 runs. 

dence on b. Figure H21 shows this dependence 
for two different values of I. For a small value 
of I, pi n u rises proportionally to b. Increas- 
ing I reveals, that, for a certain value of b, 
Pinu reaches a maximum. The value of b cor- 
responds to the slope in the log-log plot of 
the size distribution. For a larger b the prob- 
ability to find a small particle will increase. 
The summed up length of the K particles in 
the system must be equal to the length / of 
the interval for all values of b. Therefore, a 
system with a larger b has more particles but 
their average size is smaller. 

As mentioned in the previous section, 
the hierarchical initialization generates small 
holes. Thus, increasing b increases p in u 
because more smaller particles, able to fit 
into these holes, become available. On the 



11 



other hand, during each initialization, holes, 
smaller than r m j„, are generated. They will 
never be filled and increasing the number of 
particles K in the system enhances the gen- 
eration of theses holes. These two effects act 
against each other and so there exists a value 
of b for which the the initialization density 
pinit will have a maximum. 

IV. CONCLUSION 

We generalized the reversible Parking Lot 
Model introduced by Krapivsky and Ben- 
Nairn from a monodisperse model to a 
polydisperse one. This 'self consistent re- 
versible polydisperse Parking Lot Model' re- 
sults in a final and stable state of very high 
density. 

The crucial part, in order to reach high 
densities, is the self consistent desorption 
probability. The use of a constant desorption 
rate results in an exclusion of larger particles. 
The system will reach a final state of unre- 
alistic low density. Thus, the small particles 
must not be allowed to block a larger interval 
and large particles already adsorbed in the in- 
terval need to stay. This is realized through a 
hierarchical initialization where each particle 
gets I trials for adsorption. 

To keep the density on a high level we pro- 
pose a self consistent desorption probability. 
A particle can only desorb if at least one hole, 



as large as this particle, exists. The density 
rises as a power law with the number of tri- 
als / for each particle during the initialization 
and is fitted well by Eq. (J2J) . The smaller the 
range R, the lower the final density. 

The most interesting result is: There ex- 
ists a value for the exponent b of the power 
law of the to distribution for which the final 
density reaches a maximum as seen in Figure 
d 

In short: A large number of trials I and a 
wide range of particles sizes R results in a 
high final density. The maximal possible den- 
sity is determined by the exponent b. 

In this paper we extended the work of 
Brilliantov et al. [kI Q by adding the 
self consistent desorption to the polydisperse 
RSA in one dimension. We could also apply 
the expression for the evolution of density, 
measured experimentally by Knight |l| and 
Nowak to our density evolution. As the 
initialization is the crucial point in our model 
it would be interesting to see the density evo- 
lution of an experiment using our hierarchical 
initialization. 

Future work should focus on the behav- 
ior of this model in two or three dimensions 
and under the influence of gravity. The main 
difficulty will be the development of a fast 
algorithm able to determine the hole distri- 
bution. 
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